clear; clc;
%% ------------------------------------------------------------------------
% This file produces the results in Aastveit,Bjørnland and Cross (2021)
%
% This code is free to use for academic purposes only, provided that the 
% paper is cited as:
%
% Aastveit, K.A., Bjørnland H.C. and Cross, J.L. (2021). 
% Inflation expectations and the pass-through of oil prices, 
% Review of Economics and Statistics, forthcoming.
%
% This code comes without technical support of any kind.  It is expected to
% reproduce the results reported in the paper. Under no circumstances will
% the authors be held responsible for any use (or misuse) of this code in
% any way.

%% ------------------------------------------------------------------------
% Start by selecting the model specification:
% spec = 1 is the main specification
% spec = 2 is the robustness with median inflation expectations
% spec = 3 is the robustness with mixture prior on the elasticity of supply

spec = 1;

% Make directory to store results
if spec == 1
resdir = 'results_main';
disp('Replicating results for main specification')
elseif spec == 2
resdir = 'results_robust_medIE';
disp('Replicating results for robustness with median inflation expectations')
elseif spec == 3
resdir = 'results_robust_mix';
disp('Replicating results for robustness with mixture prior on elasticity of supply')
end

if ~exist(resdir, 'dir')
    mkdir(resdir)    
end

% Set paths to data, functions and stored results
addpath('data')
addpath('utilities')
addpath(resdir)

warning('off')

%% ------------------------------------------------------------------------
% Select whether to get new estimates or load stored results
 DO_NEW = 1; % Indicator variable: 1 = estimate model
%DO_NEW = 0; % Indicator variable: 0 = load stored results

% Specify inputs for MCMC
nburn=1e6;       % number of burn-in draws 
ndraws=1e5;      % number of saved MH iterations after burn-in
skip=50;         % save every "skip-th" draw to reduce autocorrelation in the Markov chain
xsi=0.04;        % tuning parameter for MH step

% Other choice parameters
nlags=12;      %number of lags
hmax=18;       %impulse response horizon (= hmax-1 months)
ndet=1;        %number of deterministic variables (1: constant)
mu=1;          %weight given to first subsample used to set prior     

% Sample
sample = (1983:1/12:2019+11/12)'; % 1983M1-2019M12 

% load oil market data from BH - sample: 1958M1 to 2016M12
data = xlsread('data_BH_AER2019','data_BH_AER2019','B3:E746'); % sample: 1958M1 to 2019M12

% column 1: global crude oil production (in million barrels/day)
% column 2: OECD+6NME industrial production
% column 3: 100*change in proxy for OECD crude oil inventories as a fraction 
%           of previous period's oil production
% column 4: real refiners' acquisition cost of crude oil imports (starts in
%           1974M1, before that '1' indicates missing data)

% transformations of oil market variables
qo=lagn(100*log(data(:,1)),1);
ip=lagn(100*log(data(:,2)),1);
stocks=data(2:end,3);
rac=lagn(100*log(data(:,4)),1); % Real IRAC
yy_rac=[qo ip rac stocks];   

% Load inflation data
if spec == 1 || 3
IE_S = xlsread('MICH','FRED Graph','C12:C515'); % Next year expectations: 1978M1-2019M12 - MEAN
elseif spec == 2
IE_S = xlsread('MICH','FRED Graph','B12:B515'); % Next year expectations: 1978M1-2019M12 - MEDIAN
end
CPI = xlsread('CPIAUCSL','FRED Graph','B12:B887'); % US CPI: 1947M1-2019M12

% Transformations of inflation data
Inf = 1200*diff(log(CPI)); % Inflation: 1947M2-2019M12

% Dates
dates_BHdat = (1958+1/12:1/12:2019+11/12)'; % Lose one obs in transformation: 1958M2 to 2019M12
dates_IE_S  = (1978:1/12:2019+11/12)'; % No transformation: 1978M1-2019M12
dates_CPI   = (1947+1/12:1/12:2019+11/12)'; % Lose one obs in transformation: 1947M2-2019M12

% Cut loaded oil market series to sample 
BH_mat = [dates_BHdat yy_rac];
ind1 = find(BH_mat(:,1) == sample(1));
ind2 = find(BH_mat(:,1) == sample(end)); 
BHdat = BH_mat(ind1:ind2,2:end); 
        
% Cut loaded inflation expectations series to sample 
IE_S_mat = [dates_IE_S IE_S];
ind1 = find(IE_S_mat(:,1) == sample(1)); 
ind2 = find(IE_S_mat(:,1) == sample(end));
InfExp = IE_S_mat(ind1:ind2,2);   
        
% Cut loaded inflation series to sample 
INF_mat = [dates_CPI Inf];
ind1 = find(INF_mat(:,1) == sample(1)); 
ind2 = find(INF_mat(:,1) == sample(end));
Inf = INF_mat(ind1:ind2,2);   

% Pool data
yy=[BHdat,InfExp,Inf];
yymat = [sample,yy]; % data with dates

%% ------------------------------------------------------------------------
% Set up model for estimation
time=sample;
n=size(yy,2);          %number of endogenous variables

%split sample
yy1=yy(1:60,:);           %first subsample 
yy2=yy(61:end,:);         %second subsample 
                           
% Get data matrices for first subsample
[X1,y1,T1]=getXy(yy1,ndet,nlags);
yyy1=y1;
xxx1=X1;

% Get data matrices for second subsample
[X2,y2,T2]=getXy(yy2,ndet,nlags);
yyy2=y2;
xxx2=X2;

seednumber=140778;
rand('seed',seednumber);
randn('seed',seednumber);           

%% ------------------------------------------------------------------------
% ALGORITHM FOR GETTING POSTERIORS OF A, B and D    
% Adapted from replication codes provided by
% Baumeister, C., & Hamilton, J. D. (2019). 
% Structural interpretation of vector autoregressions with incomplete 
% identification: Revisiting the role of oil supply and demand shocks. 
% American Economic Review, 109(5), 1873-1910.

%% ------------------------------------------------------------------------
% STEP 1a: Set parameters of the prior distributions for impact coefficients (A)

bounds = 5;
x1=-bounds:.001:0;       %grid for negative parameters
y1=0:.001:bounds;        %grid for positive parameters
z1=-bounds:.001:bounds;  %grid for parameters where no sign is imposed a priori
f1=0:.001:1;             %fraction (for beta distribution)

% alpha(qp): short-run price elasticity of oil supply (sign: positive)
if spec == 1 || spec == 2
% Student's t-distributed prior 
c_alpha_qp = 0.1;            
sigma_alpha_qp = 0.2;       
nu_alpha_qp = 3;
prior_alpha_qp=student_pos_prior(y1,c_alpha_qp,sigma_alpha_qp,nu_alpha_qp);
elseif spec == 3
% mixture prior
del4 = 0.0056;  % Taken from BH4 codes
b4 = 0.5552;    % Taken from BH4 codes
y1a=0:del4:b4;  % Grid for positive parameters
w=0.8;          % Weight on mixed prior for oil supply elasticity
c_alpha_qp = 0.1;            
sigma_alpha_qp = 0.2;       
nu_alpha_qp = 3;
es_ub = 0.0258;
prior_alpha_qp=student_pos_prior_mix_mod(y1a,c_alpha_qp,sigma_alpha_qp,nu_alpha_qp,w,es_ub);
end

% alpha(yp): short-run oil price elasticity of global demand (sign: negative)
c_alpha_yp = -0.05;
sigma_alpha_yp = 0.1;   
nu_alpha_yp = 3;
prior_alpha_yp = student_neg_prior(x1,c_alpha_yp,sigma_alpha_yp,nu_alpha_yp);

% beta(qy): income elasticity of oil demand (sign: positive)
c_beta_qy = 0.7;   
sigma_beta_qy = 0.2;   
nu_beta_qy = 3;
prior_beta_qy = student_pos_prior(y1,c_beta_qy,sigma_beta_qy,nu_beta_qy);

% beta(qp): short-run price elasticity of oil demand (sign: negative)
c_beta_qp = -0.1;   
sigma_beta_qp = 0.2;     
nu_beta_qp = 3;
prior_beta_qp = student_neg_prior(x1,c_beta_qp,sigma_beta_qp,nu_beta_qp);

% chi: OECD fraction of true oil inventories (about 60-65%)
alpha_k = 15;
beta_k = 10;
prior_k = beta_prior(f1,alpha_k,beta_k);

% psi1: short-run production elasticity of inventory demand (sign: unrestricted)
c_psi1 = 0;
sigma_psi1 = 0.5;
nu_psi1 = 3;
prior_psi1 = student_prior(z1,c_psi1,sigma_psi1,nu_psi1);

% psi3: short-run price elasticity of inventory demand (sign: unrestricted)
c_psi3 = 0;
sigma_psi3 = 0.5;
nu_psi3 = 3;
prior_psi3 = student_prior(z1,c_psi3,sigma_psi3,nu_psi3);

% rho: importance of measurement error (between 0 and chi)
% prior conditional on chi
chi = alpha_k/(alpha_k+beta_k);
mean_rho = 0.25*chi;
std_rho = 0.12*chi;
[alpha_l,beta_l] = GetBetaParameters(mean_rho,std_rho);
prior_lambda = beta_prior(f1,alpha_l,beta_l);

%Prior for det(Atilde)
zeta_h1 = 1;       % overall weight put on the prior for h1
mu_h1 = 0.6;       % location parameter of the prior (obtained by simulation: see simulate_h1.m)
sig_h1 = 1.6;      % scale parameter of the prior (obtained by simulation: see simulate_h1.m)
nu_h1 = 3;         % degrees of freedom
lam_h1 = 2;        % determines asymmetry of the prior
prior_h1 = log(1/sig_h1) + log(tpdf((z1-mu_h1)/sig_h1,nu_h1)) + log(normcdf(lam_h1*z1/sig_h1));

%Prior for H(2,2)
mu_h2 = 0.8;       % loction parameter of the prior
sig_h2 = 0.2;      % scale parameter of the prior
nu_h2 = 3;         % degrees of freedom
prior_h22 = student_prior(z1,mu_h2,sig_h2,nu_h2);

% Priors for inflation expectations
% lambda(q): oil production effect on inflation expectations (sign: negative)
c_lam_q = -.1;     
sigma_lam_q = 10;   
nu_lam_q = 3;
prior_lam_q = student_neg_prior(x1,c_lam_q,sigma_lam_q,nu_lam_q);

% gamma(y): REA effect on inflation expectations (sign: positive)
c_lam_y = 0.1; 
sigma_lam_y = 10;   
nu_lam_y = 3;
prior_lam_y = student_pos_prior(y1,c_lam_y,sigma_lam_y,nu_lam_y);

% lambda(p): oil price effect on inflation expectations (sign: positive)
c_lam_p = .02;   
sigma_lam_p = 1;   
nu_lam_p = 3;
prior_lam_p = student_pos_prior(y1,c_lam_p,sigma_lam_p,nu_lam_p);

% gamma(p): inventory effect on inflation expectations (sign: unrestricted)
c_lam_inv = 0;     
sigma_lam_inv = 100;   
nu_lam_inv = 3;
prior_lam_inv = student_prior(z1,c_lam_inv,sigma_lam_inv,nu_lam_inv);

% lambda(pi): inflation effect on inflation expectations (sign: positive)
c_lam_pi = .55;     
sigma_lam_pi = 1;   
nu_lam_pi = 3;
prior_lam_pi = student_pos_prior(y1,c_lam_pi,sigma_lam_pi,nu_lam_pi);

% Priors for inflation
% gamma(q): oil production effect on inflation (sign: negative)
c_gam_q = -.1;     
sigma_gam_q = 10;   
nu_gam_q = 3;
prior_gam_q = student_neg_prior(x1,c_gam_q,sigma_gam_q,nu_gam_q);

% gamma(y): REA effect on inflation (sign: positive)
c_gam_y = .25; % Table 3 of CGK (2018, p.1464)     
sigma_gam_y = 1;    
nu_gam_y = 3;
prior_gam_y = student_pos_prior(y1,c_gam_y,sigma_gam_y,nu_gam_y);

% gamma(p): oil price effect on inflation (sign: positive)
c_gam_p = 0.04;     
sigma_gam_p = 1;   
nu_gam_p = 3;
prior_gam_p = student_pos_prior(y1,c_gam_p,sigma_gam_p,nu_gam_p);

% gamma(p): inventory effect on inflation (sign: unrestricted)
c_gam_inv = 0;     
sigma_gam_inv = 100;   
nu_gam_inv = 3;
prior_gam_inv = student_prior(z1,c_gam_inv,sigma_gam_inv,nu_gam_inv);

% gamma(p): inflation expectations on inflation (sign: positive)
c_gam_i = 1;     
sigma_gam_i = 1;   
nu_gam_i = 3;
prior_gam_i = student_pos_prior(y1,c_gam_i,sigma_gam_i,nu_gam_i);

% Set arbitrary initial values for elements of A (prior mode/mean of
% elements in A and L)
A_old=[c_alpha_qp; c_alpha_yp; c_beta_qy; c_beta_qp; alpha_k/(alpha_k+beta_k); ...
       c_psi1; c_psi3; alpha_k/(alpha_k+beta_k)*alpha_l/(alpha_l+beta_l); ...
       c_lam_p; c_lam_pi; c_lam_y; c_lam_q; c_gam_p; c_gam_i; c_gam_y; c_gam_q];   
       
c=size(A_old,1);           %number of parameters in A to be estimated

%% ------------------------------------------------------------------------
% STEP 1b: Set informative priors on lagged coefficients (B) 
%          and for inverse of diagonal elements (D)          

% Compute standard deviation of each series residual via an OLS regression
% to be used in setting the prior (here: AR(12))
[s11,uhat1]=sd_prior(yy1(:,1),nlags);
[s22,uhat2]=sd_prior(yy1(:,2),nlags);
[s33,uhat3]=sd_prior(yy1(:,3),nlags);
[s44,uhat4]=sd_prior(yy1(:,4),nlags);
[s55,uhat5]=sd_prior(yy1(:,5),nlags);
[s66,uhat6]=sd_prior(yy1(:,6),nlags);

% See Doan (2013) for choices of values of hyperparameters (Minnesota-type prior)
lambda0=0.5;  % overall confidence in prior
lambda1=1;    % confidence on higher-order lags (lambda1 = 0 gives all lags equal weight)
lambda2=1;    % confidence in other-than-own lags 
lambda3=100;  % tightness of constant term
lambda4=100;  % tightness on lags from inflation block to oil block

% The mean for D is calibrated on diagonal elements in omega
uhat=[uhat1 uhat2 uhat3 uhat4 uhat5 uhat6];
S=uhat'*uhat/T1;

kappa=2;                          %prior 
kappastar=kappa+(mu*T1+T2)/2;     %posterior kappa 

% Specify the prior mean of the coefficients of the 6 equations of the VAR
% and their prior covariance
% PRIOR MEAN
m=zeros(n,(n*nlags+1));
m(1,3)=0.1;
m(3,3)=-0.1;

% PRIOR COVARIANCE
v1 = 1:nlags;
v1 = v1'.^(-2*lambda1);
v2 = 1./diag(S);
v3 = kron(v1,v2);
v3 = lambda0^2*[v3; lambda3^2];
v3 = 1./v3;
Pinv = diag(sqrt(v3));

% PRIOR COVARIANCE for inflation feedback to oil variables
Pinvtemp = diag(Pinv);
Pinvtemp(5:n:end-2*nlags) = lambda4; 
Pinvtemp(6:n:end-2*nlags) = lambda4; 
Pinvoil = diag(Pinvtemp);

%% ------------------------------------------------------------------------
% Complete model set up
% Compute summary statistics of the observed data for first subsample:
Syy1=yyy1'*yyy1;
Sxx1=xxx1'*xxx1;
Sxy1=yyy1'*xxx1;
zeta1=(Syy1-(Sxy1/Sxx1)*Sxy1');

% Compute summary statistics of the observed data for second subsample:
Syy2=yyy2'*yyy2;
Sxx2=xxx2'*xxx2;
Sxy2=yyy2'*xxx2;
zeta2=(Syy2-(Sxy2/Sxx2)*Sxy2');

omega_tildeT=(mu*T1+T2)\(mu*zeta1+zeta2);

Xtildeoil=[sqrt(mu).*xxx1;xxx2;Pinvoil];
Xtilde=[sqrt(mu).*xxx1;xxx2;Pinv];

% Compute M_star
M_staroil=inv(Xtildeoil'*Xtildeoil);
M_star=inv(Xtilde'*Xtilde);

%% ------------------------------------------------------------------------
% Estimation starts here
if DO_NEW
disp('Estimating Model')

% Get starting values for A via optimization 
%fixed parameter values
param=[c_alpha_qp;sigma_alpha_qp;nu_alpha_qp;c_alpha_yp;sigma_alpha_yp;nu_alpha_yp; ...
       c_beta_qy;sigma_beta_qy;nu_beta_qy;c_beta_qp;sigma_beta_qp;nu_beta_qp; ...
       alpha_k;beta_k;c_psi1;sigma_psi1;nu_psi1;c_psi3;sigma_psi3;nu_psi3; ...
       c_lam_p; sigma_lam_p; nu_lam_p; ...
       c_lam_pi; sigma_lam_pi; nu_lam_pi; ...
       c_lam_y; sigma_lam_y; nu_lam_y;...
       c_lam_q; sigma_lam_q; nu_lam_q;...
       c_gam_p; sigma_gam_p; nu_gam_p; ...
       c_gam_i; sigma_gam_i; nu_gam_i; ...
       c_gam_y; sigma_gam_y; nu_gam_y;...
       c_gam_q; sigma_gam_q; nu_gam_q;...
       vec(omega_tildeT);kappastar;kappa;T1;T2;mu_h1;sig_h1;nu_h1;lam_h1;zeta_h1; ...
       mu_h2;sig_h2;nu_h2];
   
f_anon = @(theta_hat) post_val_m_mod(theta_hat,param,S,m,yyy1,yyy2,Pinv,Pinvoil,Xtilde,Xtildeoil,mu);

%starting value for theta
theta_zero=A_old;   % mode/mean of prior distributions
options = optimset('LargeScale','off','MaxFunEvals',10000);
[theta_max,val_max,exitm,~,~,HM] = fminunc(f_anon,theta_zero,options);

% find Hessian of log posterior
if min(eig(inv(HM))) > 0
     PH = chol(inv(HM))';
else
     PH = eye(c);
end

% start MH algorithm with theta_max if sign restrictions are satisfied
if  theta_max(1,1)>0 && theta_max(2,1)<0 && theta_max(3,1)>0 && theta_max(4,1)<0 && ...
            theta_max(5,1)>0 && theta_max(5,1)<1 && theta_max(8,1)>0 && theta_max(8,1)<1 && ...
            theta_max(9,1)>0 && theta_max(10,1)>0 && theta_max(11,1)>0 && theta_max(12,1)>0 && ...
            theta_max(13,1)>0 && theta_max(14,1)>0 && theta_max(15,1)>0 && theta_max(16,1)>0
A_old=theta_max;
end

A_tilde=[1 0 -A_old(1,1) 0  0 0; ...
         0 1 -A_old(2,1) 0  0 0; ...
         1 -A_old(3:4,1)' -1/A_old(5,1)  0 0; ...
         -A_old(6,1) 0 -A_old(7,1) 1  0 0;...
         -A_old(12,1) -A_old(11,1) -A_old(9,1) 0 1 -A_old(10,1); ...
         -A_old(16,1) -A_old(15,1) -A_old(13,1) 0 -A_old(14,1) 1];
P=eye(n);
P(4,3)=A_old(8,1);
A=P*A_tilde;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%     ALGORITHM FOR GETTING POSTERIORS OF A, B and D    %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Ytilde1=[sqrt(mu).*(yyy1*A(1,:)');yyy2*A(1,:)';Pinvoil'*m(1,:)'];
Ytilde2=[sqrt(mu).*(yyy1*A(2,:)');yyy2*A(2,:)';Pinvoil'*m(2,:)'];
Ytilde3=[sqrt(mu).*(yyy1*A(3,:)');yyy2*A(3,:)';Pinvoil'*m(3,:)'];
Ytilde4=[sqrt(mu).*(yyy1*A(4,:)');yyy2*A(4,:)';Pinvoil'*m(4,:)'];
Ytilde5=[sqrt(mu).*(yyy1*A(5,:)');yyy2*A(5,:)';Pinv'*m(5,:)'];
Ytilde6=[sqrt(mu).*(yyy1*A(6,:)');yyy2*A(6,:)';Pinv'*m(6,:)'];

% Compute m_star(i)
m_star1=(Xtildeoil'*Xtildeoil)\(Xtildeoil'*Ytilde1);
m_star2=(Xtildeoil'*Xtildeoil)\(Xtildeoil'*Ytilde2);
m_star3=(Xtildeoil'*Xtildeoil)\(Xtildeoil'*Ytilde3);
m_star4=(Xtildeoil'*Xtildeoil)\(Xtildeoil'*Ytilde4);
m_star5=(Xtilde'*Xtilde)\(Xtilde'*Ytilde5);
m_star6=(Xtilde'*Xtilde)\(Xtilde'*Ytilde6);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% STEP 2: Set the variance of the candidate generating density (P) %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
W=xsi*eye(c);     %variance of RW-MH  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% STEP 3: Evaluate posterior at starting value for A:  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
omega=A*S*A';
taustar1=gettau(kappa,omega(1,1),Ytilde1,Xtildeoil);
taustar2=gettau(kappa,omega(2,2),Ytilde2,Xtildeoil);
taustar3=gettau(kappa,omega(3,3),Ytilde3,Xtildeoil);
taustar4=gettau(kappa,omega(4,4),Ytilde4,Xtildeoil);
taustar5=gettau(kappa,omega(5,5),Ytilde5,Xtilde);
taustar6=gettau(kappa,omega(6,6),Ytilde6,Xtilde);

% Evaluate prior p(A) at old draw
prior1 = student_pos_prior(A_old(1,1),c_alpha_qp,sigma_alpha_qp,nu_alpha_qp); % positive
prior2 = student_neg_prior(A_old(2,1),c_alpha_yp,sigma_alpha_yp,nu_alpha_yp); % negative
prior3 = student_pos_prior(A_old(3,1),c_beta_qy,sigma_beta_qy,nu_beta_qy); % positive
prior4 = student_neg_prior(A_old(4,1),c_beta_qp,sigma_beta_qp,nu_beta_qp); % negative
prior5 = beta_prior(A_old(5,1),alpha_k,beta_k); 
prior6 = student_prior(A_old(6,1),c_psi1,sigma_psi1,nu_psi1); % unrestricted
prior7 = student_prior(A_old(7,1),c_psi3,sigma_psi3,nu_psi3); % unrestricted
prior8 = beta_prior(A_old(8,1),alpha_l,beta_l);

prior9 = student_neg_prior(A_old(12,1),c_lam_q,sigma_lam_q,nu_lam_q); % negative
prior10 = student_pos_prior(A_old(11,1),c_lam_y,sigma_lam_y,nu_lam_y); % positive
prior11 = student_pos_prior(A_old(9,1),c_lam_p,sigma_lam_p,nu_lam_p);  % positive
prior12 = student_pos_prior(A_old(10,1),c_lam_pi,sigma_lam_pi,nu_lam_pi); % positive
prior13 = student_neg_prior(A_old(16,1),c_gam_q,sigma_gam_q,nu_gam_q); % negative
prior14 = student_pos_prior(A_old(15,1),c_gam_y,sigma_gam_y,nu_gam_y); % positive
prior15 = student_pos_prior(A_old(13,1),c_gam_p,sigma_gam_p,nu_gam_p); % positive
prior16 = student_pos_prior(A_old(14,1),c_gam_i,sigma_gam_i,nu_gam_i); % positive


h1 = det(A_tilde);
fh1 = zeta_h1*(log(tpdf((h1-mu_h1)/sig_h1,nu_h1)) + log(normcdf(lam_h1*h1/sig_h1)));

H = inv(A_tilde);
h2 = H(2,2);
prior_h2 = student_prior(h2,mu_h2,sig_h2,nu_h2);
fh2 = log(prior_h2);

% Compute posterior value at candidate draw
log_priors=log(prior1)+log(prior2)+log(prior3)+log(prior4)+ ...
           log(prior5)+log(prior6)+log(prior7)+log(prior8)+fh1+fh2+...
           +log(prior9)+log(prior10)+log(prior11)+log(prior12)+log(prior13)+log(prior14)+log(prior15)+log(prior16);
       
new_term1=kappa*log(kappa*omega(1,1));
new_term2=kappa*log(kappa*omega(2,2));
new_term3=kappa*log(kappa*omega(3,3));
new_term4=kappa*log(kappa*omega(4,4));       
new_term5=kappa*log(kappa*omega(5,5));    
new_term6=kappa*log(kappa*omega(6,6));    

new_term=new_term1+new_term2+new_term3+new_term4+new_term5+new_term6; 
up=log_priors+(mu*T1+T2)/2*log(det(A*omega_tildeT*A'))+new_term;
down=kappastar*log(2*taustar1/(mu*T1+T2)) ...
     +kappastar*log(2*taustar2/(mu*T1+T2)) ...
     +kappastar*log(2*taustar3/(mu*T1+T2)) ...
     +kappastar*log(2*taustar4/(mu*T1+T2)) ...
     +kappastar*log(2*taustar5/(mu*T1+T2)) ...
     +kappastar*log(2*taustar6/(mu*T1+T2));
posteriorOLD=up-down;

% RW-MH algorithm 
naccept=0;
count=0;

% Store posterior distributions (after burn-in)
A_post_m=zeros(c+2,ndraws);
A_post=zeros(n,n,ndraws);
B_post=zeros(n,n*nlags+1,ndraws);
D_post=zeros(n,n,ndraws);

% Misc
nstat = 0;
ndyn = 0;

while count<ndraws*skip + nburn
      count=count+1;
      if (count/10000) == floor(count/10000)
          count
      end
      
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
    % STEP 4a: Generate draw for A from the RW candidate density %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    A_new=A_old+chol(W)'*PH*randn(c,1)/sqrt(0.5*(randn(1)^2 + randn(1)^2));    % fat tails
   
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % STEP 4b: Evaluate posterior at new draw %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % impose signs/ranges
    if A_new(1,1)>0 && A_new(2,1)<0 && A_new(3,1)>0 && A_new(4,1)<0 && ...
            A_new(5,1)>0 && A_new(5,1)<1 && A_new(8,1)>0 && A_new(8,1)<1 && ...
            A_new(9,1)>0 && A_new(10,1)>0 && A_new(11,1)>0 && A_new(12,1)<0 && ...
            A_new(13,1)>0 && A_new(14,1)>0 && A_new(15,1)>0 && A_new(16,1)< 0       
        
       % Evaluate prior p(A) at new draw
        prior1 = student_pos_prior(A_new(1,1),c_alpha_qp,sigma_alpha_qp,nu_alpha_qp);
        prior2 = student_neg_prior(A_new(2,1),c_alpha_yp,sigma_alpha_yp,nu_alpha_yp);
        prior3 = student_pos_prior(A_new(3,1),c_beta_qy,sigma_beta_qy,nu_beta_qy);
        prior4 = student_neg_prior(A_new(4,1),c_beta_qp,sigma_beta_qp,nu_beta_qp);
        prior5 = beta_prior(A_new(5,1),alpha_k,beta_k);
        prior6 = student_prior(A_new(6,1),c_psi1,sigma_psi1,nu_psi1);
        prior7 = student_prior(A_new(7,1),c_psi3,sigma_psi3,nu_psi3);
        prior8 = beta_prior(A_new(8,1),alpha_l,beta_l);
        prior9 = student_neg_prior(A_new(12,1),c_lam_q,sigma_lam_q,nu_lam_q); 
        prior10 = student_pos_prior(A_new(11,1),c_lam_y,sigma_lam_y,nu_lam_y);
        prior11 = student_pos_prior(A_new(9,1),c_lam_p,sigma_lam_p,nu_lam_p); 
        prior12 = student_pos_prior(A_new(10,1),c_lam_pi,sigma_lam_pi,nu_lam_pi); 
        prior13 = student_neg_prior(A_new(16,1),c_gam_q,sigma_gam_q,nu_gam_q); 
        prior14 = student_pos_prior(A_new(15,1),c_gam_y,sigma_gam_y,nu_gam_y);
        prior15 = student_pos_prior(A_new(13,1),c_gam_p,sigma_gam_p,nu_gam_p); 
        prior16 = student_pos_prior(A_new(14,1),c_gam_i,sigma_gam_i,nu_gam_i); 

       % Construct full matrix A 
        A_tilde=[1 0 -A_new(1,1) 0  0 0; ...
                 0 1 -A_new(2,1) 0  0 0; ...
                 1 -A_new(3:4,1)' -1/A_new(5,1)  0 0; ...
                 -A_new(6,1) 0 -A_new(7,1) 1  0 0;...
                 -A_new(12,1) -A_new(11,1) -A_new(9,1) 0 1 -A_new(10,1); ...
                 -A_new(16,1) -A_new(15,1) -A_new(13,1) 0 -A_new(14,1) 1];
                 
       h1 = det(A_tilde);
       fh1 = zeta_h1*(log(tpdf((h1-mu_h1)/sig_h1,nu_h1)) + log(normcdf(lam_h1*h1/sig_h1)));

       H = inv(A_tilde);
       h2 = H(2,2);
       prior_h2 = student_prior(h2,mu_h2,sig_h2,nu_h2);
       fh2 = log(prior_h2);
            
       P=eye(n);
       P(4,3)=A_new(8,1);
       A=P*A_tilde;
       
       Ytilde1=[sqrt(mu).*(yyy1*A(1,:)');yyy2*A(1,:)';Pinvoil'*m(1,:)'];
       Ytilde2=[sqrt(mu).*(yyy1*A(2,:)');yyy2*A(2,:)';Pinvoil'*m(2,:)'];
       Ytilde3=[sqrt(mu).*(yyy1*A(3,:)');yyy2*A(3,:)';Pinvoil'*m(3,:)'];
       Ytilde4=[sqrt(mu).*(yyy1*A(4,:)');yyy2*A(4,:)';Pinvoil'*m(4,:)'];
       Ytilde5=[sqrt(mu).*(yyy1*A(5,:)');yyy2*A(5,:)';Pinv'*m(5,:)'];
       Ytilde6=[sqrt(mu).*(yyy1*A(6,:)');yyy2*A(6,:)';Pinv'*m(6,:)'];       
       
       m_star1=(Xtildeoil'*Xtildeoil)\(Xtildeoil'*Ytilde1);
       m_star2=(Xtildeoil'*Xtildeoil)\(Xtildeoil'*Ytilde2);
       m_star3=(Xtildeoil'*Xtildeoil)\(Xtildeoil'*Ytilde3);
       m_star4=(Xtildeoil'*Xtildeoil)\(Xtildeoil'*Ytilde4);
       m_star5=(Xtilde'*Xtilde)\(Xtilde'*Ytilde5);
       m_star6=(Xtilde'*Xtilde)\(Xtilde'*Ytilde6);
       
       omega=A*S*A';
       taustar1=gettau(kappa,omega(1,1),Ytilde1,Xtildeoil);
       taustar2=gettau(kappa,omega(2,2),Ytilde2,Xtildeoil);
       taustar3=gettau(kappa,omega(3,3),Ytilde3,Xtildeoil);
       taustar4=gettau(kappa,omega(4,4),Ytilde4,Xtildeoil);
       taustar5=gettau(kappa,omega(5,5),Ytilde5,Xtilde);
       taustar6=gettau(kappa,omega(6,6),Ytilde6,Xtilde);
       
       % Compute posterior value at new candidate draw
        log_priors=log(prior1)+log(prior2)+log(prior3)+log(prior4)+ ...
                   log(prior5)+log(prior6)+log(prior7)+log(prior8)+fh1+fh2+...
                   log(prior9)+log(prior10)+log(prior11)+log(prior12)+log(prior13)+log(prior14)+log(prior15)+log(prior16);
              
       new_term1=kappa*log(kappa*omega(1,1));
       new_term2=kappa*log(kappa*omega(2,2));
       new_term3=kappa*log(kappa*omega(3,3));
       new_term4=kappa*log(kappa*omega(4,4));       
       new_term5=kappa*log(kappa*omega(5,5));       
       new_term6=kappa*log(kappa*omega(6,6));       

       new_term=new_term1+new_term2+new_term3+new_term4+new_term5+new_term6;       
       up=log_priors+(mu*T1+T2)/2*log(det(A*omega_tildeT*A'))+new_term;
       down=kappastar*log(2*taustar1/(mu*T1+T2)) ...
            +kappastar*log(2*taustar2/(mu*T1+T2)) ...
            +kappastar*log(2*taustar3/(mu*T1+T2)) ...
            +kappastar*log(2*taustar4/(mu*T1+T2)) ...
            +kappastar*log(2*taustar5/(mu*T1+T2)) ...
            +kappastar*log(2*taustar6/(mu*T1+T2));
       posteriorNEW=up-down;
         
       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       % STEP 5: Compute acceptance probability %
       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       accept=min([exp(posteriorNEW-posteriorOLD);1]);
       u=rand(1);                     %draw from a uniform distribution
       if u<=accept
          A_old=A_new;                %we retain the new draw
          posteriorOLD=posteriorNEW;
          naccept=naccept+1;          %count the number of acceptances
       end
    end
    
    if count>nburn && mod( (count - nburn ), skip ) == 0
        csave = (count - nburn)/skip; % counter for storing draws

         %Store results after burn-in    
            AA_tilde=[1 0 -A_old(1,1) 0  0 0; ...
                     0 1 -A_old(2,1) 0  0 0; ...
                     1 -A_old(3:4,1)' -1/A_old(5,1)  0 0; ...
                     -A_old(6,1) 0 -A_old(7,1) 1  0 0;...
                     -A_old(12,1) -A_old(11,1) -A_old(9,1) 0 1 -A_old(10,1); ...
                     -A_old(16,1) -A_old(15,1) -A_old(13,1) 0 -A_old(14,1) 1];
             
         HH=inv(AA_tilde);               
         P=eye(n);
         P(4,3)=A_old(8,1);
         AA=P*AA_tilde;
        
         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
         % STEP 7: Generate a draw for d(ii)^-1 from independent gamma %
         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       Ytilde1=[sqrt(mu).*(yyy1*A(1,:)');yyy2*A(1,:)';Pinvoil'*m(1,:)'];
       Ytilde2=[sqrt(mu).*(yyy1*A(2,:)');yyy2*A(2,:)';Pinvoil'*m(2,:)'];
       Ytilde3=[sqrt(mu).*(yyy1*A(3,:)');yyy2*A(3,:)';Pinvoil'*m(3,:)'];
       Ytilde4=[sqrt(mu).*(yyy1*A(4,:)');yyy2*A(4,:)';Pinvoil'*m(4,:)'];
       Ytilde5=[sqrt(mu).*(yyy1*A(5,:)');yyy2*A(5,:)';Pinv'*m(5,:)'];
       Ytilde6=[sqrt(mu).*(yyy1*A(6,:)');yyy2*A(6,:)';Pinv'*m(6,:)'];       
       
       m_star1=(Xtildeoil'*Xtildeoil)\(Xtildeoil'*Ytilde1);
       m_star2=(Xtildeoil'*Xtildeoil)\(Xtildeoil'*Ytilde2);
       m_star3=(Xtildeoil'*Xtildeoil)\(Xtildeoil'*Ytilde3);
       m_star4=(Xtildeoil'*Xtildeoil)\(Xtildeoil'*Ytilde4);
       m_star5=(Xtilde'*Xtilde)\(Xtilde'*Ytilde5);
       m_star6=(Xtilde'*Xtilde)\(Xtilde'*Ytilde6);
       
       omega=AA*S*AA';
       chk = -1;
       cnt = 0;
      while chk < 0 && cnt<1000  
         cnt = cnt+1;
         d11=inv(gamrnd(kappastar,1/gettau(kappa,omega(1,1),Ytilde1,Xtildeoil)));
         d22=inv(gamrnd(kappastar,1/gettau(kappa,omega(2,2),Ytilde2,Xtildeoil)));
         d33=inv(gamrnd(kappastar,1/gettau(kappa,omega(3,3),Ytilde3,Xtildeoil)));
         d44=inv(gamrnd(kappastar,1/gettau(kappa,omega(4,4),Ytilde4,Xtildeoil)));          
         d55=inv(gamrnd(kappastar,1/gettau(kappa,omega(5,5),Ytilde5,Xtilde)));          
         d66=inv(gamrnd(kappastar,1/gettau(kappa,omega(6,6),Ytilde6,Xtilde)));          
         DD=diag([d11;d22;d33;d44;d55;d66]);
                 
         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
         % STEP 8: Generate a draw for b(i) from multivariate normal %
         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
         b1=m_star1+(randn(1,nlags*n+1)*chol(d11.*M_staroil))';
         b2=m_star2+(randn(1,nlags*n+1)*chol(d22.*M_staroil))';
         b3=m_star3+(randn(1,nlags*n+1)*chol(d33.*M_staroil))';
         b4=m_star4+(randn(1,nlags*n+1)*chol(d44.*M_staroil))';
         b5=m_star5+(randn(1,nlags*n+1)*chol(d55.*M_star))';
         b6=m_star6+(randn(1,nlags*n+1)*chol(d66.*M_star))';
         BB=[b1';b2';b3';b4';b5';b6']; 
         compBB =  varcompanion(BB,ndet,n,nlags);
            if max(abs(eig(compBB))) < 1
                chk = 1;
            end
      end
         if chk == 1
             nstat=nstat+1;
             A_post_m(:,csave)=[A_old; det(AA_tilde); HH(2,2)];    
             A_post(:,:,csave)=AA;
             D_post(:,:,csave)=DD;
             B_post(:,:,csave)=BB;
        else
             A_post_m(:,csave)=A_post_m(:,csave-1);    
             A_post(:,:,csave)=A_post(:,:,csave-1);
             D_post(:,:,csave)=D_post(:,:,csave-1);
             B_post(:,:,csave)=B_post(:,:,csave-1);            
        end
         clear AA_tilde AA BB DD       
     end 
end

% Compute acceptance ratio of RW-MH algorithm
acceptance_ratio=naccept/ndraws;
stationary_draws=nstat/ndraws;
disp(['Acceptance ratio:' num2str(acceptance_ratio)])
disp(['Stationary draws ratio:' num2str(stationary_draws)])

save(sprintf('%s/results_posterior_draws.mat',resdir), 'A_post_m', 'A_post', 'D_post', 'B_post', '-v7.3')

else
disp('Loading stored draws from joint posterior')
    load(sprintf('%s/results_posterior_draws.mat',resdir))    
end

%% ------------------------------------------------------------------------
% Generate Figures
if spec == 1
% plot the prior and posterior distributions for the elements in A and H
figure_posteriors

% Analyze convergence of the chain
p1=0.1;   %first 10% of the sample (for Geweke's (1992) convergence diagnostic)
p2=0.5;   %last 50% of the sample
names = {'${\alpha}_{qp}$',...
         '${\alpha}_{yp}$',...
         '${\beta}_{qy}$',...
         '${\beta}_{qp}$',...
         '${\chi}$',...
         '${\psi}_{1}$',...
         '${\psi}_{3}$',...
         '${\rho}$',...
         '${\lambda}_{\pi^e p}$',...
         '${\lambda}_{\pi^e \pi}$',...
         '${\lambda}_{\pi^e y}$',...
         '${\lambda}_{\pi^e q}$',...
         '${\gamma}_{\pi p}$',...
         '${\gamma}_{\pi \pi^e}$',...
         '${\gamma}_{\pi y}$',...
         '${\gamma}_{\pi q}$',...
         'det(${\tilde A}$)',...
         '${h}_{22}$'};
disp('Plotting Figure A7: Traceplots of MCMC draws, and generating Table A3: Chi-squared probabilities for equality of means using Geweke Chi-squared test for each parameter chain')    
[autoc,Table_chi2] = convergence_diagnostics(A_post_m,p1,p2,names,resdir);
save(sprintf('%s/results_conv_diags.mat',resdir), 'autoc', 'Table_chi2')
end

figure_HDECs_compute % computes the historical decompositions
figure_IRFs_compute  % computes actual and counterfactual IRFs

figure_IRFs_plot_inf % plot actual IRFs for inflation block  
figure_IRFs_plot_CF  % plot counterfactual IRFs for inflation block
figure_GR_CF         % plot Great recession counterfactual

if spec == 1
figure_HDECs_plot    % plots the historical decompositions
figure_IRFs_plot_oil % plot actual IRFs for oil block       
KStests     % Two-sample Kolmogorov-Smirnov goodness-of-fit hypothesis test 
            % for difference between actual and counterfactual IRFs
end

